################################################
#### Read in the Odds ratios and their CIS. ####

library(convertGraph)

matched_or <- read.csv("/Users/sharpe/Desktop/Young Surgeons/taper_table/binary_outcomes_1.2.2019.csv")
level_3_matched_or <- read.csv("/Users/sharpe/Desktop/Young Surgeons/table/bin_outcomes_era_by_surg_REDONE_11.26.2018.csv")
level_3_matched_or_both <- read.csv("/Users/sharpe/Desktop/Young Surgeons/table/bin_outcomes_trad_mod_REDONE_11.26.2018.csv")
level_3_matched_or <- subset(level_3_matched_or, comparison %in% c("mod gensurg","mod ortho"))
level_3_matched_or_both <- subset(level_3_matched_or_both, comparison=="mod")

both_labels <- c("Level 0 - Gensurg + Ortho", "Level 1 - Gensurg + Ortho", "Level 2 - Gensurg + Ortho")
gensurg_labels <- c("Level 0 - Gensurg", "Level 1 - Gensurg", "Level 2 - Gensurg")
ortho_labels <- c("Level 0 - Ortho", "Level 1 - Ortho", "Level 2 - Ortho")

zero_labels <- c("Level 0 - Gensurg + Ortho", "Level 0 - Gensurg", "Level 0 - Ortho")
one_labels <- c("Level 1 - Gensurg + Ortho", "Level 1 - Gensurg", "Level 1 - Ortho")
two_labels <- c("Level 2 - Gensurg + Ortho", "Level 2 - Gensurg", "Level 2 - Ortho")

matched_or$surg_type <- ifelse(matched_or$comparison %in% both_labels, "Both", ifelse(matched_or$comparison %in% gensurg_labels, "Gensurg", "Ortho"))
matched_or$level <- ifelse(matched_or$comparison %in% zero_labels, 0, ifelse(matched_or$comparison %in% one_labels, 1, ifelse(matched_or$comparison %in% two_labels, 2, 3)))

level_3_matched_or$surg_type <- ifelse(grepl(pattern="gensurg", x=level_3_matched_or$comparison), "Gensurg", "Ortho")
level_3_matched_or$level <- 3

level_3_matched_or_both$surg_type <- "Both"
level_3_matched_or_both$level <- 3 

matched_or <- subset(matched_or, select=c("Outcome","comparison","surg_type", "level", "Paired.Odds.Ratio",  "Odds.Ratio.Lower.95.", "Odds.Ratio.Upper.95."))
level_3_matched_or <- subset(level_3_matched_or, select=c("Outcome","comparison","surg_type", "level", "Paired.Odds.Ratio",  "Odds.Ratio.Lower.95.", "Odds.Ratio.Upper.95."))
level_3_matched_or_both <- subset(level_3_matched_or_both, select=c("Outcome","comparison","surg_type", "level", "Paired.Odds.Ratio",  "Odds.Ratio.Lower.95.", "Odds.Ratio.Upper.95."))

all_data <- rbind(matched_or, level_3_matched_or, level_3_matched_or_both)
all_data <- all_data[order(all_data$surg_type, all_data$level, decreasing=T),]

##################################
#### Read in the M-estimates. ####

matched_m <- read.csv("/Users/sharpe/Desktop/Young Surgeons/taper_table/matched_m_estimates_1.2.2019_FIXED.csv")
matched_m_lv_3 <- read.csv("/Users/sharpe/Desktop/Young Surgeons/table/simul_matched_m_estimates_REDONE_11.28.2018.csv")

matched_m$surg_type <- ifelse(matched_m$Comparison %in% both_labels, "Both", ifelse(matched_m$Comparison %in% gensurg_labels, "Gensurg", "Ortho"))
matched_m$level <- ifelse(matched_m$Comparison %in% zero_labels, 0, ifelse(matched_m$Comparison %in% one_labels, 1, ifelse(matched_m$Comparison %in% two_labels, 2, 3)))

matched_m_lv_3 <- subset(matched_m_lv_3, Comparison %in% c("mod", "mod gensurg", "mod ortho"))
matched_m_lv_3$level <- 3
matched_m_lv_3$surg_type <- ifelse(matched_m_lv_3$Comparison=="mod", "Both", ifelse(matched_m_lv_3$Comparison=="mod gensurg", "Gensurg", "Ortho"))

m_data <- rbind(matched_m, matched_m_lv_3)
m_data <- m_data[order(m_data$surg_type, m_data$level, decreasing=T),]

##############################################################################################
#### Create Figure 1, which consists of mortality, readmission, plos, anesthia time, FTR. ####

# pdf(file="/Users/sharpe/Desktop/Young Surgeons/taper_table/dead_30_OR_plot.pdf", width=8.5)
# pdf(file="/Users/sharpe/Desktop/Young Surgeons/taper_table/figure_1.pdf", width=8.5)
tiff(file="/Users/sharpe/Desktop/Young Surgeons/taper_table/figure_1_v1_FIXED_4.1.2019.tiff", height=7, width=8, units='in', res=300)



vars_to_plot <- c("dead_30", "ftr_30_use_poa", "readms_or_death_new", "prolonged_d")
vars_to_plot <- rev(vars_to_plot)
paper_names <- c("30-day all          \n location Mortality", "30-day                 \n Failure-to-Rescue","30-day                       \n Readmission or Death","Prolonged Length of Stay")
paper_names <- rev(paper_names)
# var_names <- c("30-day Mortality", "Readmission or ", "")

all_data_2 <- subset(all_data, Outcome %in% vars_to_plot)
all_data_2$y_axis <- ifelse(all_data_2$Outcome==vars_to_plot[1], 1,
                            ifelse(all_data_2$Outcome==vars_to_plot[2], 2,
                                   ifelse(all_data_2$Outcome==vars_to_plot[3], 3, 
                                          ifelse(all_data_2$Outcome==vars_to_plot[4], 4,
                                                 ifelse(all_data_2$Outcome==vars_to_plot[5], 5, 6)))))

graph_colors <- c("gray70","gray50","gray30","gray0")
stagger <- c(0.08,0.23,0.38,0.53)
graph_labels <- c("Baseline","Baseline/Operations","Baseline/Operations/Emergent","Baseline/Operaitons/Emergent/Risk")

layout(matrix(c(1,1,1,1,
                1,1,1,1,
                1,1,1,1,
                1,1,1,1,
                1,1,1,1,
                2,2,2,2,
                2,2,2,2), nrow=7, byrow=T))

for(i in 1:length(vars_to_plot)){
  
  temp_data <- subset(all_data_2, surg_type=="Both" & Outcome==vars_to_plot[i])
  
  if(i==1){
    par(mar=c(5,5,0.2,2))
    plot(x=all_data_2$Paired.Odds.Ratio, y=all_data_2$y_axis, xlim=c(0.8,1.6), xaxt="n", yaxt="n", ylab="", xlab="Paired Odds Ratios for New Surgeon versus Experienced Surgeon", type="n", ylim=c(1,4.8))
    axis(side=1, at=c(seq(0.8,1.6,0.1)), labels=c(seq(0.8,1.6,0.1)))
    abline(h=c(1.9,2.9,3.9), lty=2)
    abline(v=c(1), lty=2)
    text(x=0.85, y=1.40, labels=paper_names[i], cex=1)
    text(x=0.84, y=2.45, labels=paper_names[i+1], cex=1)
    text(x=0.83, y=3.45, labels=paper_names[i+2], cex=1)
    text(x=0.83, y=4.4, labels=paper_names[i+3], cex=1)
    # legend(x=1.45, y=4, col=rev(graph_colors), legend=c("Level 3", "Level 2", "Level 1", "Level 0"), lty=c(1,1,1,1))
  }
  
  for(j in 3:0){
    print(j)
    
    temp_temp_data <- subset(temp_data, level==j)
    temp_temp_data$y_axis <- temp_temp_data$y_axis+stagger[j+1]
    points(x=temp_temp_data$Paired.Odds.Ratio, y=temp_temp_data$y_axis, col=graph_colors[j+1], lty=1, pch=19)
    lines(x=c(temp_temp_data$Odds.Ratio.Lower.95., temp_temp_data$Odds.Ratio.Upper.95.), y=c(i+stagger[j+1],i+stagger[j+1]), col=graph_colors[j+1])
    lines(x=c(temp_temp_data$Odds.Ratio.Lower.95., temp_temp_data$Odds.Ratio.Lower.95.), y=c(i+stagger[j+1]-0.05, i+stagger[j+1]+0.05), col=graph_colors[j+1])
    lines(x=c(temp_temp_data$Odds.Ratio.Upper.95., temp_temp_data$Odds.Ratio.Upper.95.), y=c(i+stagger[j+1]-0.05, i+stagger[j+1]+0.05), col=graph_colors[j+1])
    
  }
  
  
}



vars_to_plot_m <- c("anesth_time_claim")
paper_names_m <- c("Anesthesia Time (minutes)")
m_data_2 <- subset(m_data, Outcome %in% vars_to_plot_m)
m_data_2$y_axis <- 1

graph_colors <- c("gray70","gray50","gray30","gray0")
stagger <- c(0.12,0.28,0.42,0.57)


for(i in 1:length(vars_to_plot_m)){
  
  temp_data <- subset(m_data_2, surg_type=="Both" & Outcome==vars_to_plot_m[i])
  
  if(i==1){
    par(mar=c(8,5,0,2))
    plot(x=m_data_2$New...Exp..M.estimate, y=m_data_2$y_axis, xlim=c(-6.65,20), xaxt="n", yaxt="n", ylab="", xlab="Paired Difference (M-Estimate) for New Surgeon versus Experienced Surgeon", type="n", ylim=c(1.00,1.9))
    axis(side=1, at=c(seq(-5,20,1)), labels=c(seq(-5,20,1)))
    text(x=-4.6, y=1.45, labels=paper_names_m[i], cex=1)
    abline(v=c(0), lty=3)
    par(xpd=T)
    legend(x=3, y=0.4, col=rev(graph_colors), legend=c("Baseline + Operations + Emergent + Risk Factors", "Baseline + Operations + Emergent", "Baseline + Operations", "Baseline"), lty=c(1,1,1,1), horiz=F, box.lty=0,
           cex=0.85)
  }
  
  for(j in 3:0){
    print(j)
    
    temp_temp_data <- subset(temp_data, level==j)
    temp_temp_data$y_axis <- temp_temp_data$y_axis+stagger[j+1]
    points(x=temp_temp_data$New...Exp..M.estimate, y=temp_temp_data$y_axis, col=graph_colors[j+1], lty=1, pch=19)
    lines(x=c(temp_temp_data$New...Exp..Lower.95..CI, temp_temp_data$New...Exp..Upper.95..CI), y=c(i+stagger[j+1],i+stagger[j+1]), col=graph_colors[j+1])
    lines(x=c(temp_temp_data$New...Exp..Lower.95..CI, temp_temp_data$New...Exp..Lower.95..CI), y=c(i+stagger[j+1]-0.035, i+stagger[j+1]+0.035), col=graph_colors[j+1])
    lines(x=c(temp_temp_data$New...Exp..Upper.95..CI, temp_temp_data$New...Exp..Upper.95..CI), y=c(i+stagger[j+1]-0.035, i+stagger[j+1]+0.035), col=graph_colors[j+1])
    
  }
  
  
}

dev.off()

#############################################################################################
#### Create figure 2, which contains ICU Usage, Length of Stay and 30-day Resource Cost. ####

# pdf(file="/Users/sharpe/Desktop/Young Surgeons/taper_table/figure_2.pdf", width=8.5)
tiff(file="/Users/sharpe/Desktop/Young Surgeons/taper_table/figure_2_v1_FIXED_4.1.2019.tiff", height=7, width=8, units='in', res=300)

vars_to_plot <- c("icu_yes")
paper_names <- c("ICU Usage")

all_data_3 <- subset(all_data, Outcome %in% vars_to_plot)
all_data_3$y_axis <- 1

graph_colors <- c("gray70","gray50","gray30","gray0")
stagger <- c(0.17,0.32,0.47,0.62)

layout(matrix(c(1,1,1,1,
                1,1,1,1,
                2,2,2,2,
                2,2,2,2,
                3,3,3,3,
                3,3,3,3), nrow=6, byrow=T))

for(i in 1:length(vars_to_plot)){
  
  temp_data <- subset(all_data_3, surg_type=="Both" & Outcome==vars_to_plot[i])
  
  if(i==1){
    par(mar=c(6,5,0.2,2))
    plot(x=all_data_3$Paired.Odds.Ratio, y=all_data_3$y_axis, xlim=c(0.8,1.6), xaxt="n", yaxt="n", ylab="", xlab="Paired Odds Ratios for New Surgeon versus Experienced Surgeon", type="n", ylim=c(1,1.8))
    axis(side=1, at=c(seq(0.8,1.6,0.05)), labels=c(seq(0.8,1.6,0.05)))
    text(x=0.81, y=1.40, labels=paper_names[i], cex=1)
    abline(v=c(1), lty=3)
  }
  
  for(j in 3:0){
    print(j)
    
    temp_temp_data <- subset(temp_data, level==j)
    temp_temp_data$y_axis <- temp_temp_data$y_axis+stagger[j+1]
    points(x=temp_temp_data$Paired.Odds.Ratio, y=temp_temp_data$y_axis, col=graph_colors[j+1], lty=1, pch=19)
    lines(x=c(temp_temp_data$Odds.Ratio.Lower.95., temp_temp_data$Odds.Ratio.Upper.95.), y=c(i+stagger[j+1],i+stagger[j+1]), col=graph_colors[j+1])
    lines(x=c(temp_temp_data$Odds.Ratio.Lower.95., temp_temp_data$Odds.Ratio.Lower.95.), y=c(i+stagger[j+1]-0.035, i+stagger[j+1]+0.035), col=graph_colors[j+1])
    lines(x=c(temp_temp_data$Odds.Ratio.Upper.95., temp_temp_data$Odds.Ratio.Upper.95.), y=c(i+stagger[j+1]-0.035, i+stagger[j+1]+0.035), col=graph_colors[j+1])
    
  }
  
  
}




vars_to_plot_m <- c("los_long")
paper_names_m <- c("Length of Stay (Days)")
m_data_3 <- subset(m_data, Outcome %in% vars_to_plot_m)
m_data_3$y_axis <- 1

graph_colors <- c("gray70","gray50","gray30","gray0")
stagger <- c(0.08,0.23,0.38,0.53)


for(i in 1:length(vars_to_plot_m)){
  
  temp_data <- subset(m_data_3, surg_type=="Both" & Outcome==vars_to_plot_m[i])
  
  if(i==1){
    par(mar=c(6,5,0.2,2))
    plot(x=m_data_3$New...Exp..M.estimate, y=m_data_3$y_axis, xlim=c(-0.366,1.1), xaxt="n", yaxt="n", ylab="", xlab="Paired Difference (M-Estimate) for New Surgeon versus Experienced Surgeon", type="n", ylim=c(0.98,1.7))
    axis(side=1, at=c(seq(-0.3,-0.1,0.1),0,seq(0.1,1.1,0.1)), labels=c(seq(-0.3,-0.1,0.1),0,seq(0.1,1.1,0.1)))
    text(x=-0.29, y=1.35, labels=paper_names_m[i], cex=1)
    abline(v=0, lty=3)
    # par(xpd=T)
  }
  
  for(j in 3:0){
    print(j)
    
    temp_temp_data <- subset(temp_data, level==j)
    temp_temp_data$y_axis <- temp_temp_data$y_axis+stagger[j+1]
    points(x=temp_temp_data$New...Exp..M.estimate, y=temp_temp_data$y_axis, col=graph_colors[j+1], lty=1, pch=19)
    lines(x=c(temp_temp_data$New...Exp..Lower.95..CI, temp_temp_data$New...Exp..Upper.95..CI), y=c(i+stagger[j+1],i+stagger[j+1]), col=graph_colors[j+1])
    lines(x=c(temp_temp_data$New...Exp..Lower.95..CI, temp_temp_data$New...Exp..Lower.95..CI), y=c(i+stagger[j+1]-0.03, i+stagger[j+1]+0.03), col=graph_colors[j+1])
    lines(x=c(temp_temp_data$New...Exp..Upper.95..CI, temp_temp_data$New...Exp..Upper.95..CI), y=c(i+stagger[j+1]-0.03, i+stagger[j+1]+0.03), col=graph_colors[j+1])
    
  }
  
  
}



vars_to_plot_m <- c("total_cost_complete")
paper_names_m <- c("30-day Resource Cost\n(2013 US Dollars)       ")
m_data_3 <- subset(m_data, Outcome %in% vars_to_plot_m)
m_data_3$y_axis <- 1

graph_colors <- c("gray70","gray50","gray30","gray0")
stagger <- c(0.09,0.24,0.39,0.52)


for(i in 1:length(vars_to_plot_m)){
  
  temp_data <- subset(m_data_3, surg_type=="Both" & Outcome==vars_to_plot_m[i])
  
  if(i==1){
    par(mar=c(8,5,0.0,2))
    plot(x=m_data_3$New...Exp..M.estimate, y=m_data_3$y_axis, xlim=c(-1065,3200), xaxt="n", yaxt="n", ylab="", xlab="Paired Difference (M-Estimate) for New Surgeon versus Experienced Surgeon", type="n", ylim=c(0.98,1.7))
    axis(side=1, at=c(seq(-1100,-200,300),0,seq(200,3200,300)), labels=c(seq(-1100,-200,300),0,seq(200,3200,300)))
    text(x=-790, y=1.35, labels=paper_names_m[i], cex=1)
    abline(v=c(0), lty=3)
    par(xpd=T)
    legend(x=450, y=0.62, col=rev(graph_colors), legend=c("Baseline + Operations + Emergent + Risk Factors", "Baseline + Operations + Emergent", "Baseline + Operations", "Baseline"), lty=c(1,1,1,1), horiz=F, box.lty=0,
    cex=0.85)
}
  
  for(j in 3:0){
    print(j)
    
    temp_temp_data <- subset(temp_data, level==j)
    temp_temp_data$y_axis <- temp_temp_data$y_axis+stagger[j+1]
    points(x=temp_temp_data$New...Exp..M.estimate, y=temp_temp_data$y_axis, col=graph_colors[j+1], lty=1, pch=19)
    lines(x=c(temp_temp_data$New...Exp..Lower.95..CI, temp_temp_data$New...Exp..Upper.95..CI), y=c(i+stagger[j+1],i+stagger[j+1]), col=graph_colors[j+1])
    lines(x=c(temp_temp_data$New...Exp..Lower.95..CI, temp_temp_data$New...Exp..Lower.95..CI), y=c(i+stagger[j+1]-0.03, i+stagger[j+1]+0.03), col=graph_colors[j+1])
    lines(x=c(temp_temp_data$New...Exp..Upper.95..CI, temp_temp_data$New...Exp..Upper.95..CI), y=c(i+stagger[j+1]-0.03, i+stagger[j+1]+0.03), col=graph_colors[j+1])
    
  }
  
  
}


dev.off()



